function [regions]=assembleHistoneDSforClustering()
% Organizes histone modification data (GM06990, chrm 7) by chromosome regions
%
% assembleHistoneDSforClustering() constructs a Region view of the data, by
% reorganizing histone modification data saved on 'GM06990chips_chrm7.mat'.
% The new DS 'regions' contains histone measurements organized by regions
% of the chromosome 7, and is used by performHistoneClustering() as data 
% source to plot the cluster tree and heatmap.
% 
% 'regions' also integrates data from different DAS sources. Annotation of
% genes is fetched from the Havana Gene Database at Sanger Institute, while
% the cancer link is obtained from the Genetic Association Database (GAD)
% at NCBI. 
%
% NOTE: Assembling regions can take several hours, since thousands of URL calls
% are made to retrieve information at Havana and GAD databases.
%
% SYNTAX: [] = assembleHistoneDS()
%
% See also assembleHistoneDS(), performHistoneClustering_Fig4CD()
%  
%    DASMiner: DAS library and browser for Matlab.
%    Diogo Veiga, Jan 2009.

load('GM06990chips_chrm7.mat', 'histoneChips');

regions = struct('label',{},'samples',{},'values',{} );
regCount = 0;


for i=1:size(histoneChips.chips,2)
   
    for j=1:size(histoneChips.chips(i).region,2)
       
        regLabel = ['chr' histoneChips.chips(i).region(j).chrm ':' histoneChips.chips(i).region(j).start ...
            ',' histoneChips.chips(i).region(j).stop];
        
        res = strcmpi(regLabel,{regions.label});
        if (any(res))
            index = find(res,1); %get index for label matching
            
            regions(index).samples = [regions(index).samples ; histoneChips.chips(i).sample];
            regions(index).values = [regions(index).values ; histoneChips.chips(i).values(j)];
        
        else %add a new region
           regCount = regCount+1; 
           regions(regCount).label = regLabel;
           regions(regCount).samples = {histoneChips.chips(i).sample};
           regions(regCount).values = histoneChips.chips(i).values(j);
           regions(regCount).geneSymbol = annotateGeneName(histoneChips.chips(i).region(j).chrm, ...
               histoneChips.chips(i).region(j).start, histoneChips.chips(i).region(j).stop);
        end
        
    end
end

for i=1:size(regions,2)
    regions(i).cancerTypes = annotateCancerLink(regions(i).geneSymbol);
end

save GM06990chips_chrm7_RegionsView.mat regions;

function [geneSymbol] = annotateGeneName(chrm,start,stop)
%Accessing Havana Genes Database
%eg. http://das.sanger.ac.uk/das/otter_das/features?segment=7:116391273,116392269

[xmlS] = executeDASCommand('http://das.sanger.ac.uk/das','command', 'features', 'DSN', 'otter_das','chrom',chrm,'start',start ...
    , 'stop', stop, 'timeout', 5000);

geneSymbol = '';
if (isstruct(xmlS)) % ow problem in the server 
    if ( isfield(xmlS.GFF{1}.SEGMENT{1},'FEATURE')    )

        for i=1:size(xmlS.GFF{1}.SEGMENT{1}.FEATURE,2)
            if (~strcmp (xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.NOTE{7}.CONTENT(10:end),geneSymbol))
                if (strcmp(geneSymbol,''))
                    geneSymbol = [xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.NOTE{7}.CONTENT(10:end) ','];
                elseif (strfind(geneSymbol,[xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.NOTE{7}.CONTENT(10:end) ','])) %already on the list
                    continue;
                else
                    geneSymbol = [geneSymbol xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.NOTE{7}.CONTENT(10:end) ',']; %new gene
                end
            end
        end

    end
end

function [cancerTypes] = annotateCancerLink(geneSymbol)

cancerTypes = '';
geneTokens = regexp(geneSymbol,',','split');

for k=1:(size(geneTokens,2)-1)
    % Accessing the Genetic Association Database
    % eg. http://www.ebi.ac.uk/das-srv/genedas/das/gad/features?segment=GSTT1
    [xmlS] = executeDASCommand('http://www.ebi.ac.uk/das-srv/genedas/das','command', 'features', 'DSN', 'gad','segment',geneTokens{k});

 if (isstruct(xmlS)) % ow problem in the server  
    if ( isfield(xmlS.GFF{1}.SEGMENT{1},'FEATURE') )

        for i=1:size(xmlS.GFF{1}.SEGMENT{1}.FEATURE,2)

            for j=1:size(xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.NOTE,2)

                if (strfind(xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.NOTE{j}.CONTENT,'CANCER/'))
%                     cancerTypes = [cancerTypes ',' xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.NOTE{j}.CONTENT(8:end)];
                    if (strcmp(cancerTypes,''))
                        cancerTypes = [xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.NOTE{j}.CONTENT(8:end) ','];
                    elseif (strfind(cancerTypes,[xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.NOTE{j}.CONTENT(8:end) ','])) %already on the list
                        continue;
                    else
                        cancerTypes = [cancerTypes xmlS.GFF{1}.SEGMENT{1}.FEATURE{i}.NOTE{j}.CONTENT(8:end) ',']; %new gene
                    end
                end
            end
            
        end
    end
 end

end